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Abstract —Contingency screening for transient stability of large 
scale, strongly nonlinear, interconnected power systems is one of 
the most computationally challenging parts of Dynamic Security 
Assessment and requires huge resources to perform time-domain 
simulations-based assessment. To reduce computational cost of 
time-domain simulations, direct energy methods have been exten¬ 
sively developed. However, these methods, as well as other existing 
methods, still rely on time-consuming numerical integration of 
the fault-on dynamics. This task is computationally hard, since 
possibly thousands of contingencies need to be scanned and 
thousands of accompanied fault-on dynamics simulations need 
to be performed and stored on a regular basis. In this paper, we 
introduce a novel framework to eliminate the need for fault-on 
dynamics simulations in contingency screening. This simulation- 
free framework is based on bounding the fault-on dynamics and 
extending the recently introduced Lyapunov Function Family 
approach for transient stability analysis of structure-preserving 
model. In turn, a lower bound of the critical clearing time (CCT) 
is obtained by solving convex optimization problems without re¬ 
lying on any time-domain simulations. A comprehensive analysis 
is carried out to validate this novel technique on a number of 
IEEE test cases. 

I. Introduction 

Transient stability assessment, concerned with power sys¬ 
tems stability/instability after contingencies, is a core element 
of the Dynamic Security Assessment Systems monitoring and 
allowing the reliable operation of power systems around the 
world. The most straightforward and dominant approach in 
industry to this problem is based on the direct time-domain 
simulations of transient post-fault dynamics following possible 
contingencies. Rapid advances in computational hardware en¬ 
able it to perform accurate simulations of large scale systems 
possibly faster than real-time 0, El- However, in practice 
there are usually thousands to millions of contingencies that 
need to be screened on a regular basis. As such, the com¬ 
putational cost for time-domain simulations-based transient 
stability assessment is huge. At the same time, most of these 
contingencies are not critical, and thus most of computational 
resources are spent for assessment of contingencies that do 
not contribute to overall system risk. 

To avoid time-consuming numerical integration of post-fault 
dynamics and save the computational resources, the smarter 

Thanh Long Vu and Konstantin Turitsyn are with the Department of 
Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, 
MA, 02139 USA, e-mail: longvu@mit.edu and turitsyn@mit.edu. Surour A1 
Araifi and Mohamed Elmoursi are with Department of Electrical Engineer¬ 
ing and Computer Science, Masdar Institute, Abu Dhabi, U.A.E., email: 
salaraifi@masdar.ac.ae and melmoursi@masdar.ac.ae. 


way nowadays is to use a combination of the direct energy 
approaches and time-domain simulation 0-El, in which 
most contingencies will be screened by the energy method 
and the remaining contingencies are checked by time-domain 
simulations. The advantage of direct energy method is that it 
allows fast screening of contingencies while providing math¬ 
ematically rigorous certificates of stability. After decades of 
research and development, the controlling unstable equilibrium 
point (UEP) method 0 has been widely accepted as the 
most successful method among other energy function based 
direct screening methods, and is being applied in industry. 
This method is based on comparing the post-fault energy with 
the energy at the controlling UEP to certify transient stability. 

The noticeable drawback of the controlling UEP method is 
the inherent difficulty of directly identifying the controlling 
UEP 0. The controlling UEP is defined as the first UEP 
whose stable manifold is hit by the fault-on trajectory at the 
exit point, i.e. the point where the fault-on trajectory meets the 
actual stability boundary of the post-fault Stable Equilibrium 
Point (SEP). Note that the actual stability boundary of the 
SEP is generally unknown, and thus the computation of the 
exit point is very complicated and usually necessitates iterative 
time-domain simulations. For a given fault-on trajectory, the 
controlling UEP computation requires solving a large set of 
nonlinear differential algebraic equations which is done by 
numerical methods. However, with respect to these methods, 
e.g. Newton method, the convergence region of the controlling 
UEP can be very small and irregular compared to that of the 
SEP. If an initial guess for the numerical solver was not suf¬ 
ficiently close to the controlling UEP, then the computational 
algorithm will result in wrong controlling UEP and might 
probably converge to a SEP, leading to unreliable stability 
assessment. Unfortunately, it is extremely hard to find an initial 
guess sufficiently close to the controlling UEP. 

The second drawback of the controlling UEP method is that 
it requires simulating and storing each fault-on trajectory to 
carry out the assessment for the respective contingencies. To 
the best of our knowledge, there are only a few works on 
contingency screening without relying on fault-on dynamics 
simulations. Particularly, in 0 the closest UEP method is 
exploited and an algebraic formulation of the critical clearing 
time is obtained based on polynomial approximation of the 
swing equations. However it is assumed that the dynamics 
of the rotor angles during the fault is a constant positive 
acceleration. This approximation is remarkable and may cause 
incorrect estimation of the critical clearing time. 

The objective of this paper is to develop novel numeri- 
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cal approach that can potentially alleviate the computational 
burden of finding the controlling UEP. We aim to achieve 
this objective by developing a completely simulation-free 
technique for the estimation of critical clearing time. This 
technique is based on an extension of the recently introduced 
Lyapunov Functions Family (LFF) approach 0. The principle 
of this approach is to provide transient stability certificates by 
constructing a family of Lyapunov functions and then finding 
the best suited function in the family for given initial states. 
Basically, this method certifies that the post-fault dynamics 
is stable if the fault-cleared state stays within a polytope 
surrounding the post-fault equilibrium point and the Lyapunov 
function at the fault-cleared state is smaller than the minimum 
value of Lyapunov function over the flow-out boundary of that 
polytope. Therefore, to screen the contingencies for transient 
stability, this method only requires the knowledge of the fault- 
cleared state, instead of the whole fault-on trajectory. 

Exploiting this advantage of LFF method, a technique is 
introduced to bound the fault-on dynamics and thereby the 
fault-cleared state. This bound leads to a transient stability 
certificate that only relies on checking the clearing time, i.e. 
if the clearing time is under certain threshold then the fault- 
cleared state is still in the region of attraction of the original 
SEP and the post-fault dynamics is determined stable. By this 
new method, a fast transient stability assessment for a large 
number of contingencies can be obtained without using any 
simulations. Such approach can be utilized in several power 
system applications, such as optimal power flow, resources 
allocation, and HVDC control problems ma-im where 
the proposed transient stability certificate can help reduce 
the search space by eliminating less critical contingencies in 
studies. 

The structure of this paper is as follows. In Section [TT| 
the contingency screening problem addressed in this paper is 
introduced, together with the extension of the LFF approach 
for transient stability analysis. Section [III] presents the main 
result of this paper regarding the simulation-free algebraic 
estimation of the critical clearing time, and explains how this 
new stability certificate can be used in practice to screen 
contingency for transient stability without any time-domain 
simulations. Finally, in Section [TV] performance of the pro¬ 
posed method on contingency screening of several IEEE test 
systems is presented and analyzed. Section [V] concludes the 
paper with discussions about possible ways to improve the 
algorithms. 

II. Lyapunov Function Family Approach for 
Transient Stability 

In this section, we show that the Lyapunov function family 
approach 0, originally presented for the Kron-reduction 
model, is applicable to the transient stability analysis of 
structure-preserving power models. Then, we extend this fam¬ 
ily to a set of convex Lyapunov functions family, that will 
be instrumental to establish a lower bound of critical clearing 
time in the next section. 

In normal conditions, power grids operate at some stable 
equilibrium point. During disturbances such as faults, the 


system evolves subject to the fault-on (disturbance) dynamics 
and moves away from the pre-fault equilibrium point. After 
the fault is cleared, the system may return back to the pre¬ 
fault SEP or to a new post-fault SEP depending on whether the 
fault is self-cleared or cleared by circuit breakers action. In this 
paper, the proposed method tackles the type of contingencies, 
where a fault occurs in a transmission line and then self 
clears such that the post-fault network recovers to the pre¬ 
fault network topology. To describe the post-fault dynamics, 
we utilize the differential structure-preserving model fl8l . This 
model naturally incorporates the dynamics of rotor angle as 
well as response of dynamic load power output to frequency 
deviation. Though it does not model the dynamics of voltage 
in the system, in comparison to the Kron-reduction models 
with constant impedance loads ED, the structure of power 
systems and the impact of load dynamics are preserved in 
this approach. When the losses of the transmission lines are 
ignored, the model can be expressed as: 


m k S k + d k S k + ^2 a k? sin (4 - 4) =Pm k i (1) 

j&J^k 

k = 1 ,..., to , 

d k S k T ^ ' d'kj sin(4 Sj ) — (2) 




k = m + 1,..., n, 


where the first m equations represent the dynamics of gen¬ 
erators and the remaining (n — m) equations represent the 
dynamics of frequency-dependent loads. With k = 1 
then m k is the dimensionless moment of inertia of the k th 
generator, d k is the term representing primary frequency 
controller action on the governor, and P rrik is the effective 
dimensionless mechanical power input acting on the rotor. 
With k = m + 1,..., n , then d k > 0 is the constant frequency 
coefficient of load and P^ k is the nominal load. Let £ be the 
set of all the transmission lines and Af k be the set of neigh¬ 
boring buses of the bus k th . Then, a k j = V k VjB k j , where 
[B k j]{ k j}££ is the susceptance matrix and V k represents the 
voltage magnitude at the k th bus, both of which are assumed 
to be constant. The stationary operating condition is given by 
[4,..., <5*, 0,..., 0] T where 6 k is solution of the power flow¬ 
like equations 

VI a kj sin(4 - Sj) = P k , Vfc = 1,..., n, (3) 

where P k = P mk ,k = and P k = ~P° k ,k = 

m + 1,..., n. We assume that there exists a stable operating 
condition S* £ A(A),A < 7r/2, where the polytope A(A) is 
defined by inequalities \S k j\ < A for all {k, j} £ £. 

In the LFF approach, the nonlinear couplings and the 
linear model are separated. To do that, the state vector 
x = [xi, X 2 , X 3 ] t is introduced which is composed of 
the vector of generator’s angle deviations from equilibrium 
x\ = [4 — 4,..., 5 m — S* J T , their angular velocities 
X 2 = [4,... ,S m ] T , and the vector of load’s angle deviation 
from equilibrium X 3 = [4i+i — S ^ +1 ,..., 4 — 5^} T . Let 
E be the incidence matrix of the corresponding graph, so 
that E[Si...S n ] T = [(4 — 4){fc,j}ef] T - Consider matrix 


IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , NOV. 2015 


3 



Fig. 1. Strict bounding of the nonlinear function by linear functions of 
the angular difference .) in the lossy power systems 


C such that Cx = E[5±... 5 n ] 7 . Consider the vector of 
nonlinear power flow F in the simple trigonometric form 
F{Cx ) = [(sin 5kj - smS* kj ) {k}j}&£ ] T . 

Then, in state space representation the system can be 
expressed in the following compact form: 


Xi = x 2 

x 2 = M^D lX2 - S l D~ 1 E T SF{Cx) (4) 

±3 = -S 2 D~ 1 E T SF{Cx) 

where S = diag [o-kj){k,j}e£ is the diagonal matrix of 
coupling magnitudes and Si = [/ mxm O mxn - m ],S 2 = 

[On—mxm ^n—mxn—m\^E \ = diag(di, . . . , d m ), A'/i = 

diag(m 1; ..., m n ), D = diag(mi,... , m m , d m +i, ■ ■ ■ , d n ). 
Equivalently, 

x = Ax — BF(Cx), (5) 


with the matrices A, B given by the following expression: 



Omxm 

1 mxm 

Omxn—m 

A = 

Omxm 


Omxn—m 

and 

On—mxm 

On—mxm 

O n—mxn—m 

B 


S 1 D~ 1 E t S ; 

s 2 d~ 1 e t s ] . 


Here, |£j is the number of edges in the graph defined by the 
susceptance matrix, or equivalently the number of non-zero 
non-diagonal entries in B k j. 

For the system defined by (0, the LFF approach proposes 
to use the Fyapunov functions family given by: 


V(x) = -x t Qx — 


E 

{k,j}e£ 


K 


{k,j } (cos S kj + S k j sin 5* kj ) 


(7) 


in which the diagonal, nonnegative matrices K, H and the 
symmetric, nonnegative matrix Q satisfy the following linear 
matrix inequality (FMI): 


A t Q + QA R 
R t -2H 


( 8 ) 


with R = QB — C T H — (KCA) T . Then, it can be proved that 
the Fyapunov function is decreasing in the polytope V defined 


by inequalities \S k j+S k j\ < n,\/{k,j} € 8. In order to ensure 
that the system will not escape the polytope V during transient 
dynamics one condition will be added to restrict the set of 
initial states inside V. Accordingly, we define the minimization 
of the function V(x) over the union dV out of the flow-out 
boundary segments dV as follows: 

Kin = min V(x), (9) 

x£dV out 

where dV is the flow-out boundary segment of polytope 
V that is defined, for each transmission line {k,j} € 8 
connecting generator buses k and j, by \S k j + S* kj \ = 7r and 
^kj Skj > 0. Given the value of I4iin, an FFF-based estimation 
for the region of attraction of the equilibrium point is given 
by 

n v = {x g v : V{x) < H min }. (10) 

Finally, to determine if the post-fault dynamics is stable, we 
check if the fault-cleared state xo is inside the stability region 
estimate IZ-p, i.e. if xo is in the polytope V and V (xo) < Vmin- 
Therefore, to certify transient stability of each contingency, the 
FFF approach only need to know the fault-cleared state xq (i.e. 
the state of fault-on trajectory at the clearing time), rather than 
the whole fault-on trajectory. 

In this paper, the proposed approach is only concerned with 
voltage phase angles staying inside the polytope Q defined 
by inequalities |<5^^-1 < n/2,\/{k,j} € 8. An advantage of 
considering this polytope of voltage phasor angles is that 
inside this polytope the Fyapunov function V(x) defined in 0 
is convex. As such, the minimum value V m i n can be calculated 
in polynomial time. In addition, inside this polytope, a stricter 
bounding for the nonlinear flow vector F can be established 
as follows 


(f{kj} - {S kj - s* kj ))(f {ktj} - P(Skj - S* kj )) <0 ( 11 ) 


where j3 = 


1 — sin A 


> 0 and f{k,j} = sin Skj — sin 5 k 


is 


7t/2 — X ' - J ' KJ> . KJ 

an element of the vector F. Exploiting this strict bound of the 

nonlinear flow vector F, the FMI 0 can be replaced by the 
following less restrictive FMI: 


' A t Q + QA- 2/3C t HC R 
R t ~2H 

R = QB - (1 + (3)C t H - ( KCA) t , 


( 12 ) 


while all the above results for the stability certificate still hold 
true. In particular, the estimate for region of attraction is given 
by 

K q = {x G Q : V{x) < 1W (13) 

with 


Kin = min V{x). (14) 

:Dsae»“‘ 

The proof of this fact is given in Annendix IVI-AI With the less 
restrictive FMI dT2l> . a broader family of Fyapunov functions 
can be obtained, which will be exploited to establish the lower 
bound of the critical clearing time in the next section. 

Remark 1: The main drawback of the proposed stability 
certificate is that it currently does not incorporate the detailed 
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model of generators and its associate control systems, such as 
excitation systems, PSS and governor system. Swing equation 
model doesn’t incorporate associated control systems and 
generator’s fast dynamics and assumes a fixed field voltage 
magnitude during transient period. However, the setpoint val¬ 
ues of voltage magnitude can be allowed to fluctuate around 
the nominal value Vo (let say less than 10% around Vo). In 
the matrix B in ([ 6 ]), we take the new the coupling magnitude 
diagonal matrix S = diag(l.l 2 VfB k j)^ k jy G g. Consider the 
new nonlinear vector F = [fkj]{k,j}e£ where 


fkj 


V k Vj (sin 5 k j - sinJjy) 


(15) 


We can see that the bounding for nonlinear function f k j in 
<DD still holds true with /3 replaced by the smaller value 
0.9 2 /3/l.l 2 . Then, all the other results will follow accordingly. 
As such, the simple Lyapunov function and stability 
region estimate ( 10 can be easily extended to the case when 
voltage magnitude setpoints fluctuate 10 % around the nominal 
value. In this case, since we have looser bounding for the 
nonlinear vector F, the according stability region estimate will 
be smaller than the original defined in (10 . Therefore, the pro¬ 
posed framework can manifest the fact that the stability region 
is smaller due to the effects of generators’ control systems 
(however, it cannot capture the voltage collapse phenomenon 
when the voltage magnitudes sag to the low values). From 
this analysis, we suggest that in the practical transient stability 
assessment, we should accordingly modify the estimation of 
the stability region to avoid overestimation of the CCT due to 
the usage of simple generators’ model. 

Remark 2: Since the proposed stability certificate only 
requires the Lyapunov function to be locally decreasing, rather 
than decreasing in the whole state space as in the energy 
method, the LFF framework can be extended to incorporate the 
losses in transmission lines. Indeed, the stability analysis here 
is essentially based on bounding the nonlinear function fkj by 
linear functions of S k j as in (fill , i.e. whenever the bounding 
(DU holds true, we can have the stability region estimate ac¬ 
cordingly. For the power systems with losses, we take the cou¬ 
pling magnitude diagonal matrix S = diag (V k VjYkj) { k ,j}^£ 
and the nonlinear function fkj as 


fkj = (sin (Skj + akj) - sin($jL + a kj ) (16) 

Here, Y kj = G 2 kj + B\. and a k j = arctan (G kj /B kj ) < 1 , 

where G k j and B k j are the (normalized) conductance and 
susceptance of the transmission line {k, j}. From Fig. Q] we 
can show that the nonlinear bounding (fTTb still holds true for 

any x £ Q and 


/3 = min 


sin(7r/2 + a kj ) - sin(|<Sjy| + otkj) 


r/2- | (5; 


(17) 


kj\ 


Then, all the stability analysis follows accordingly. Therefore, 
the LFF framework and the CCT estimation to be presented 
in the next section is applicable to lossy power systems. We 
will illustrate the proposed framework for estimating CCT of 
the lossy 2-bus system in Section IV.A. 


III. Contingency Screening without Time-domain 
Simulations 

In this section, we present a new approach to the contin¬ 
gency screening problem, which relies on a combination the 
LFF framework introduced in the previous section and the 
bounding for the reachability set of the fault-on dynamics, 
through which we can guarantee that the fault-cleared state 
is still inside the region of attraction of the post-fault stable 
equilibrium point. Interestingly, this bound leads to an alge¬ 
braic simulation-free lower bound of the critical clearing time. 
Therefore, this contingency screening approach completely 
removes any time-domain simulations of both the post-fault 
dynamics and fault-on dynamics. 


A. Bounding for The Fault-on Dynamics 

If the time-domain simulation for fault-on dynamics is used, 
the fault-cleared state xo can be determined by directly inte¬ 
grating the fault-on dynamics. Then, the value of Vq = V(xq) 
computed from 0 is compared to the value of Vmin to certify 
transient stability. 

Now, assume that time-domain simulations are not used to 
integrate the fault-on dynamics. Then the fault-cleared state 
x'o will not be known precisely. To guarantee that xq £ Q and 
V (xq) < Vmin, we will bound the fault-on dynamics. Consider 
the normal condition when the pre-fault system is in the stable 
operating condition defined by the stable equilibrium point 
5* re £ A (A). Assume that a fault occurs at the transmission 
line [ u. v} £ £ and then self-clears such that the power 
network recovers to its pre-fault topology. During the fault, 
the power system dynamics is approximated by equations: 


Xp — Axp BF pre (Cxp) -t- BD^ U sin(5 U np (18) 


Here, the fault-on trajectory is denoted as Xf ( t ) to differentiate 
it from the post-fault trajectory x(t) in 0. is the 

unit vector to extract the nonlinear function (sin S UVF — 
sin S * lVpre ) from the nonlinear vector F pre = [(sin S k j F — 
sin S k j )]{fc,j}e£> which serves to model the elimination of 
the faulted line {it, u} during the fault. In Appendix IVI-BI the 
following center result regarding the bounding of the fault- 
on dynamics is proven, which will be instrumental to the 
introduction of stability certificate in the next section. If there 
exist matrices Q, K 1 H, H > 0 and a positive number 7 such 
that 

A + ^(QBD^ uv y)(QBD^ uv y) T 
R t 


R 

-2H 


< 0, (19) 


where A = A T Q+QA-2/3C T HC, R = QB- ( 1+/3)C T H - 
(I\CA ) T , then along the fault-on dynamics ( fl8] > we have 


V(xf( t)) < — whenever xp(t) being in the polytope Q. 

27 

Note that due to <®, the Lyapunov function’s derivative 
V ( x) along the post-fault dynamics 0 is non-positive in the 
polytope Q. Basically, the above result provides a certificate 
to make sure that the fault-on dynamics does not deviate too 
much from the post-fault dynamics. As such, if the clearing 
time is under some threshold, then the fault-cleared state (i.e. 
the state of fault-on system at the clearing time) is not very 
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Fig. 2. Algorithm to screen contingencies for transient stability without 
simulations of fault-on dynamics and post-fault dynamics 


far from the considered working condition. The above result 
as such is essential to establish a lower bound of the critical 
clearing time in the next section. 

B. Estimation of The Critical Clearing Time 

Let the clearing time be T c i ear ing- In Appendix IVI-C1 the 
following stability certificate which only relies on checking 
the clearing time is proven. If the inequality ( 1 1 9| ) holds and 
the clearing time T c i earing satisfies Tearing < 2')(V m ia - 
V(xp re )), where x pre = 6* re — 5 * ost , then, the fault-cleared 
state xf{T c i e aring) is still inside the region of attraction of the 
post-fault SEP S* ost and the post-fault dynamics following the 
considered contingency leads to the stable operating condition 
A* 

w post ' 

Therefore, this stability certificate provides us with a lower 
bound of the critical clearing time as 2 r y(V m i n — V(x pre )) 
obtained by solving the inequality i fThl i. This estimation is 
totally simulation-free, distinguishing it from other methods 
in the literature to estimate the critical clearing time. 

We note that it is also possible to extend this stability 
certificate to the case when several contingencies co-exist. This 


case is of practical interest. Indeed, the large-area blackout in 
practice is usually a result of multiple contingencies happening 
at short time interval. Though large-area blackout is rare, its 
effect is severe, both economically and humanly. Therefore, 
it is critical to check if the power grids stand when several 
contingencies are happening, or leading to large-area blackout. 
The technique presented in this paper provides a framework 
to certify the safety of power grids. 


C. Choosing Lyapunov Function and Parameter 7 

Since there is a family of Lyapunov functions V(x). charac¬ 
terized by matrices Q,K , and positive numbers 7 that satisfy 
the inequality (fl9l) . we have different estimations 2^(V m m — 
V ( x pre )) °f the critical clearing time (CCT). To get the highest 
possible estimation of the CCT, we need to find the maximum 
value of 27 (Vmin — V(x P re)) over all the matrices Q,K and 
positive numbers 7 satisfying C3). Unfortunately, this is an 
NP-hard, strongly nonlinear optimization problem with both 
nonlinear objective function and nonlinear constraint. 

We observe that a good selection of Lyapunov function and 
the parameter 7 is obtained if we can predict the location 
of the fault-cleared state. In the following, we propose two 
procedures suggesting some directions to search for feasible 
Lyapunov function and parameter 7 allowing for good estima¬ 
tion of the CCT. The first procedure is totally heuristic, where 
we vary 7 and find the corresponding Lyapunov function. 
The second one is based on a prediction of the fault-cleared 
state. Both of these procedures rely on solving a number of 
convex optimization problems in the form of either quadratic 
programming or semidefinite programming. 

Procedure 1: To solve the inequality (IT9l ). we note that for 
a fixed value of 7 , the inequality Cl3 can be transformed 
to the following LMI of the matrices Q,K,H via Schur 
complement: 

' A T Q + QA-2PC T HC (^ 7 (QBD {u , v} ) R) 1 
(V7 (QBD {uM ) R) t -L \ ~ ’ 

( 20 ) 


where L = 


I 

O 


O 
2 H 


. The matrices Q. K. H can be found 


quickly from the LMI (l20t by convex optimization. Therefore, 
a heuristic algorithm can be used to find solution of (IT9l >. in 
which 7 is varied and the LMI (l2fTb is solved to obtain the 
matrices Q , K , H accordingly. 

Procedure 2: 

1) Calculate the distance r from the equilibrium point S post 
to the boundary of the polytope Q as r = minaggg ||A — 

Spost II2- 

2) Take k points x\ ,...,Xk uniformly distributed on the 
sphere S = {A : ||A — < 5* ost ||2 = r} which surrounds 
S post and stays inside Q. These points are considered as 
possible predictions for the fault-cleared state. 

3) For each point X{, using the adaptation algorithm pro¬ 
posed in J9l to find a Lyapunov function I 7 (x) character¬ 
ized by matrices (}, , Ki such that the point . 7 ;,; stays inside 
the stability region estimate 7 Zq defined in m. This 
adaptation algorithm can quickly find a suitable Lyapunov 
function after a finite number of steps. 
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4) For the matrices Qi,Ki , find the maximum value 7* 
satisfying the inequality as: 7* = max 7 subject 
to (TPjl ) where Q = Qi.K = Ki,H = Hi. Calculate 

n = 2 7 *(y mini - Vi(x pre )). 

5) Take the estimation of the CCT as the maximum value 
out of Tl, Tfc. 

We note that compared to Procedure 1, Procedure 2 may 
remarkably increase the computational complexity of calcu¬ 
lating the CCT estimate. Recent studies shown that matrices 
appearing in power system context are characterized by graphs 
with low maximal clique order, and thus the related SDP in 
these procedures can be quickly solved by the new gener¬ 
ation of SDP solvers l20il . I2T1 . In addition, the advances 
in parallel computing, e.g. distributed computing with zero 
overhead communication, promises to significantly reduce the 
computational load for these SDP solvers. 

D. Contingency Screening without Simulations 

The stability certificate in Section III.B provides us with 
a way to directly screen contingencies for transient stability 
assessment without any time-domain simulations, as described 
by the algorithm in Fig. [2] Basically, for the contingency 
manifested by the tripping of line {u, v}, one can check 
if the inequality flU) is solvable. In case it is solvable to 
find the matrices Q,K, H, and the positive number 7, then 
the Lyapunov function V(x) can be derived as in ©, and 
the minimum value V m \ u defined in (IT4ll can be calculated. 
Finally, if the clearing time (CT) T c i ear i ng satisfies that 

Tclearing 4 27(I / min ^ (*Tpre))? where Xp re &pre ^posti 

then we conclude that the post-fault dynamics following the 
considered contingency leads to a stable operating condition. 
If this inequality is not true, or if there is no solution for 
the inequality (IT9l) . then nothing can be concluded about 
the stability or instability of the post-fault dynamics. The 
contingency in this case should be screened by other energy 
method or by direct time-domain simulations. 

In contingency screening, it is greatly advantageous if we 
have a certificate to screen any possible contingency associated 
with the tripping of any transmission line in the set T C £. 
Let D be a matrix larger than or equal to D{ u ,v}D^ u „} f° r 
all the lines {u,u} £ T. We have the following result for 
the robust screening of contingencies. If the inequality (fl9t 
holds with D{ u v yDj u ^ replaced by D, and the clearing 
time Tciearing satisfies T c i ear i n g <C 27 (Vj n i n V(xp re )), then, 
for any contingency associated with the tripping of any line 
{u, u} £ IF, the fault-cleared state XF{T c iearing) is still 
inside the region of attraction of the post-fault SEP 5* ost , and 
the post-fault dynamics following the considered contingency 
leads to the stable operating condition 5* ost . This result is a 
straightforward corollary of the stability certificate in Section 
IIII-B1 and thus its proof is omitted here. 

IV. Numerical Illustrations 
A. Classical 2-Bus lossy System with Different Pre-fault and 
Post-fault SEPs 

For illustrating the presented concepts, this section presents 
the simulation results on the most simple 2-bus lossy power 



Fig. 3. System trajectory according to the fault-on dynamics and post-fault 
dynamics with the clearing time CT = — V{x vre )) = 1.0600s 



Fig. 4. Dynamics of the Lyapunov function during the fault-on stage and post¬ 
fault stage with the clearing time CT = 27 (V r m i n — V(x pr e)) = 1.0600s 


system, described by the single 2-nd order differential equation 

mS + dS + asin(i5 + a) — p = 0. (21) 

For numerical simulations, we choose m = 0.1 p.u., d = 0.15 
p.u., a = 0.2 p.u., and a = 0.05 rad. The pre-fault and post¬ 
fault power inputs are p pre = 0.05 p.u. and p pos t = 0.06 
p.u. Then, the pre-fault and post-fault stable equilibrium point 
are given by [S* re 0] T = [0.2027 0] T and [S* ost 0] T = 
[0.2547 0] T , both of which are in the polytope A(7r/10). 
Hence, /3 = (sin(7r/2 + a) — sin(7r/10 + a))/(it/2 — 7r/10) = 
0.5114. By varying 7 and solving the LMI (l20b . we obtain 
the corresponding lower bounds for the critical clearing time 
as in Tab. [J 


7 

^(I'lnin VX#pre))(s) 

1 

0.9442 

2 

0.9757 

3 

1.0077 

4 

1.0297 

5 

1.0439 

6 

1.0535 

7 

1.0600 

8 

1.0578 

9 

1.0574 

10 

1.0553 


TABLE I 

Lower bound of the critical clearing time vs. 7 
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Therefore, in these values of 7, with 7 = 7 we obtain the 
largest lower bound for the critical clearing time as 1.0600. 
The corresponding matrices Q,K,H are 


Q 


0.0443 0.0127 
0.0127 0.0879 


; K = 0.0968; H = 0.2412, (22) 


while the corresponding value of Vjnin — V(x pre ) is 0.0528. 
In Fig. [3] we show the dynamics of the system trajectory in 
the fault-on and post-fault-stage in which the clearing time is 
taken as T c i earin g = 2^{V min -V (x pre )) = 1.0600s. It can be 
seen that when the fault happens, the system evolves according 
to the fault-on dynamics and the system trajectory deviates 
from the pre-fault equilibrium point 8 * re to the fault-cleared 
state 8 fault-cleared- After the fault self-clears, the system 
trajectory recovers from the fault-cleared state 8 fault-cleared, 
to the post-fault equilibrium point 8 * t which is different 
from the pre-fault equilibrium. Figure [4] shows the divergence 
of the Lyapunov function during the fault-on stage and the 
convergence of Lyapunov function during the post-fault stage. 
These figures confirm the estimation of the critical clearing 
time as obtained by the proposed method in this paper. 


B. Three Generator System 

Consider the system of three generators with the time- 
invariant terminal voltages and mechanical torques given in 
Tab. [n] 


Node 

V (p.u.) 

Pk (P-u.) 

1 

1.0566 

-0.2464 

2 

1.0502 

0.2086 

3 

1.0170 

0.0378 


TABLE II 

Voltage and mechanical input 


The susceptances of the transmission lines are B 12 = 0.739 
p.u., .B 13 = 1.0958 p.u., and B 23 = 1.245 p.u. The equilib¬ 
rium point is calculated as: <5* = [—0.6634 — 0.5046 — 
0.5640 0 0 0 ] r , which belongs to the polytope A( 7 t/ 10 ). 
Hence, we can take /3 = (1 — sin( 7 r/ 10 ))/( 7 r /2 — 7 r/ 10 ). For 
simplicity we just take rrik = 2, dk = 1, k = 1,2, 3. Assume 
that the fault happens at the transmission line connecting 
generators 1 and 2 and then self-clears. Also, during that time 
the mechanical inputs are assumed to be unchanged. Taking 
7 = 3 and using CVX software we can solve the LMI (l20l > 
we obtain Q as 


3.8376 

3.8012 

3.5779 

7.5549 

7.4619 

7.4166 

3.8012 

3.8457 

3.5698 

7.4776 

7.5530 

7.4029 

3.5779 

3.5698 

4.0690 

7.4010 

7.4185 

7.6140 

7.5549 

7.4776 

7.4010 

38.9402 

38.2449 

38.0704 

7.4619 

7.5530 

7.4185 

38.2449 

38.9534 

38.0571 

7.4166 

7.4029 

7.6140 

38.0704 

38.0571 

39.1280 


(23) 


and K = diag(0.2554, 0.3638, 0.4386), H = 

diag(0.0943, 0.2533, 0.2960). The corresponding estimation of 
the critical clearing time is — V(x pre )) = 0.2376s. 



C. Kundur 9-Bus 3-Generator System 

Consider the Kundur 9 bus 3 machine system depicted 
in Fig. 0 with 3 generator buses and 6 load buses. 
The susceptances of the transmission lines are as fol¬ 
lows: B 14 = 17.3611p.Ti., B 27 = 16.0000p.u., H 39 = 
17.0648p.7t., -B 45 = ll.7647p.Tt., H 57 = 6.2112p.Tt., Bq^ = 
10.8696p.7t., -B 78 = 13.8889p.7t., Hgg = 9.9206p.Tt., Bqq = 
5.8824p.Ti. The bus voltages 14, mechanical inputs P mk , and 
steady state load P® are given in Tab. [HI] The stable 
operating condition is obtained by solving equations ([3]) as 
x* = [0.0381 0.3208 0.1924 - 0.0349 - 0.0421 - 

0.0409 0.0519 0.0178 0.0155 00000000 0], 
which stays in the polytope A(7 t/8 ). Hence >3 = (1 - 
sin( 7 r/ 8 ))/( 7 r /2 — 7 t/8 ) = 0.5240. The parameters for gen¬ 
erators are mi = 0.1254, m 2 = 0.034, m 3 = 0.016, d\ = 
0.0627, d 2 = 0.017, ^3 = 0.008. For simplicity, we take 
dk = 0.05, k = 4..., 9. Assume that the fault trips the line 
between buses 6 and 4 and when the fault is cleared this 
line is re-closed. With 7 = 7.10 -6 , using the CVX software, 
we can solve the LMI (l20l) in Is to obtain the Lyapunov 
function. Accordingly, we can calculate the minimum value 
of the Lyapunov function and obtain the estimation for the 
critical clearing time as 27(I4iin — V(x pre )) = 0.1175s. 


Node 

V (p.u.) 

Pk (P-u.) 

1 

1.0284 

0.6700 

2 

1.0085 

1.6300 

3 

0.9522 

0.8500 

4 

1.0627 

-0.5000 

5 

1.0707 

-0.7500 

6 

1.0749 

-0.4500 

7 

1.0490 

-0.4500 

8 

1.0579 

-0.5000 

9 

1.0521 

-0.5000 


TABLE III 

Bus voltages, mechanical inputs and static loads 


We perform time domain simulations to find the critical 
clearing time for the system when the generators are modeled 
by swing equations and by 4 th orders machine models incor¬ 
porating generators’ control systems. Accordingly, we can find 
that when the fault happens at the transmission line {4, 6 }, the 
true critical clearing times for the swing model and 4 th orders 
machine models are, respectively, 0.25s and 0.18s. Therefore, 
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the critical clearing time estimated by the proposed method in 
this paper is about half of the true one. We conclude that the 
proposed method is conservative in comparison to the time 
domain simulations, but there is no overestimation for the 
CCT. In addition, the time domain simulations confirm the 
analysis we described in Remark 1 that the generators’ control 
systems make the critical clearing time to reduce. 

In comparison to the controlling UEP method, the proposed 
method in this paper is also more conservative since the 
controlling UEP was reported 0 to get the estimate for critical 
clearing time which is different in less than 10% from the 
true one obtained by time domain simulation. However, we 
note that the CCT estimate proposed in this paper does not 
require time-domain simulation for the fault-on dynamics as 
in the controlling UEP method. This will help significantly 
reduce the computational resources spent for contingency 
screening. Therefore, the proposed framework in this paper can 
be considered as a complement of the time domain simulation 
method and controlling UEP method, which could be effi¬ 
ciently used when we aim to screen non-critic al contingencies 
with little computational resources. 

V. Conclusions and Path Forward 

In this paper, we introduced techniques to screen contin¬ 
gencies for transient stability without relying on any time- 
domain simulations. This is based on extending the recently 
introduced LFF transient stability certificate in the combina¬ 
tion with bounding for the fault-on dynamics. Basically, the 
LFF approach can certify the post-fault dynamics’s stability 
when the fault-cleared state is in some polytope surrounding 
the post-fault stable operating point and the Lyapunov function 
at the fault-cleared state is under some threshold. We observed 
that the LFF certificate only needs to know the fault-cleared 
state, instead of the fault-on trajectory. Therefore, with the 
introduced bounding technique we can bound the Lyapunov 
function at the fault-cleared state, by which we certify sta¬ 
bility for a given contingency scenario without involving any 
simulations for the fault-on trajectory and post-fault trajectory. 
In turns, we obtained an algebraic formulation for the lower 
bound of the critical clearing time, and hence the stability 
assessment only involved checking if the clearing time is 
smaller than that lower bound to assure the stability of the 
post-fault dynamics. Remarkably, the proposed stability cer¬ 
tificate only relies on solving convex optimization problems. 
It may be therefore scalable to contingency screening of 
large scale power systems, especially when combined with 
the recent advances in semi-definite programming exploiting 
the relatively low tree-width of the grids’ graph l20l . 

Toward the practical applications of the proposed 
simulation-free approach to contingency screening, further 
extensions should be made in the future where more compli¬ 
cated models of power systems and faults are considered, e.g. 
generators’ control systems, effects of buses’ reactive power, 
and permanent faults are incorporated. First, since the LFF 
method is applicable to lossy power grid lf22l . it is possible 
to extend the proposed method in this paper to incorporating 
reactive power, which will introduce the cosine term in the 


model (0. This can be done by extending the state vector 
x and combining the technique in this paper with the LFF 
transient stability techniques for lossy power grids (without 
reactive power considered) l22l . Second, we can see that, in 
order to make sure the Lyapunov function is decreasing in 
the polytope Q , it is not necessary to restrict the nonlinear 
terms F(Cx ) to be univariate. As such, we can extend the 
proposed method to power systems with generators’ voltage 
dynamics in which the voltage variable is incorporated in a 
multivariable nonlinear function F. Last, the important class 
of permanent faults, which will also result in non-identical 
pre-fault and post-fault SEPs, should be considered in the 
future work with further mathematical development for the 
representation of system dynamics under faults and more 
sophisticated estimation of critical clearing time. 

In the applications, the proposed simulation-free contin¬ 
gency screening method could be developed to robustly assess 
the stability of power systems when a set of faults happen. This 
will be applicable to assess major blackout. Also, such a robust 
certificate can be applied when there are significant changes 
in the power gird topology such as in load shedding |[23ll- 
|25| and controlled islanding schemes 126l-li30|. For this end, 
a more restrictive bounding of the fault-on dynamics should 
be employed to alleviate the conservativeness of the proposed 
method which is expected when multiple faults are considered. 

VI. Appendix 

A. Proof of the Transient Stability Certificate 

From the inequality (Ell. there exist matrices 

-^|£| x(n+m )5 Y\£\ x \s\ such that 

A t Q + QA- 2f3C T HC = - X T X : 

QB - (1 + f5)C T H - (KCA) t = - X t Y 1 

-2H = - Y t Y. 

The derivative of V(x) along 0 is hence given by: 

V(x) = 0.5 x t Qx + 0.5x t Qx 

- ^2 K{ k jy(- sin S k j + sin S kj )6 k j 

= 0.5 x T (A T Q + QA)x — x t QBF + F T KCx 

= 0.5x T (2f3C T HC - X T X)x 

- x T ((l + /3)C t H + (. KCA) t - X t Y)F 

+ F T KC(Ax - BF ) (24) 

Noting that CB = 0 and Y T Y = 2H yields 

V(x) = -0.5(Xx - YFf(Xx -YF) + J2 H {k , j }g { k, j } 

(25) 

where g {kJ} = (f {ktj} - (6 kj -<%•))(/{*,j} ~P(foj ~ S kj)) < 
0,Vx € Q. As such, the Lyapunov function V(x) is decaying 
inside the polytope Q. The other results immediately follow 
those in 0. 
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B. Proof of the Bounding of Fault-on Dynamics 

From the inequality ( fl9l ), there exist matrices 

X\£\ x(n+m)j Y\g\ x \s\ such that 

A t Q + QA- 2 PC T HC + 7 {QBD UV )(QBD UV ) T = - X T X , 
QB - (1 + p)C t H - ( KCA) t = - X t Y, 

-2 H = - Y t Y. 
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